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Abstract: The dynamics of a model, originally proposed for a type of insta- 
bility in plastic flow, has been investigated in detail. The bifurcation portrait 
of the system in two physically relevant parameters exhibits a rich variety of 
dynamical behaviour, including period bubbling and period adding or Farey 
sequences. The complex bifurcation sequences, characterized by Mixed Mode 
Oscillations, exhibit partial features of Shilnikov and Gavrilov-Shilnikov sce- 
nario. Utilizing the fact that the model has disparate time scales of dynam- 
ics, we explain the origin of the relaxation oscillations using the geometrical 
structure of the bent-slow manifold. Based on a local analysis, we calculate 
the maximum number of small amplitude oscillations, s, in the periodic orbit 
of L s type, for a given value of the control parameter. This further leads to a 
scaling relation for the small amplitude oscillations. The incomplete approach 
to homoclinicity is shown to be a result of the finite rate of 'softening' of 
the eigen values of the saddle focus fixed point. The latter is a consequence 
of the physically relevant constraint of the system which translates into the 
occurrence of back-to-back Hopf bifurcation. 
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1 Introduction 



Number of autonomous dynamical systems exhibit complex bifurcation se- 
quences with alternate periodic-chaotic states in control parameter space. The 
chaotic states are usually predominant mixture of the periodic states occur- 
ring on either side of the chaotic states. These periodic states are characterized 
by combinations of relatively large amplitude excursions and small amplitude 
near harmonic oscillations of the trajectories and have been referred to as 
mixed mode oscillations (MMOs) in the literature. The MMOs and the ac- 
companying complex bifurcation sequences have been observed in models and 
experiments in various fields of chemical kinetics fl|-f|, electrochemical reac- 
tions [pHT^I , biological systems O], and in many other physical systems [|T^,p~5 



Both numerical as well as analytical studies have been carried out extensively 
to the explain the origin of the MMOs and the complex bifurcation sequences 
these systems exhibit |]9|- [T2| , p!6| -^] . Even though the origin of complex bifur- 
cation sequences and the accompanying MMOs may depend on the particular 
system under study, almost all proposed mechanisms suggest that these bi- 
furcation portraits are the artefact of global bifurcations of the system |fl7f1 . 
Investigations into the global nature of the bifurcations of the MMOs and the 
complex bifurcation sequences have shown that homoclinic bifurcations may 
be relevant to a wide variety of systems which display the MMOs. ShilnikovE| 



has shown that if a dynamical system possesses a homoclinic orbit which is 
bi-asymptotic to a saddle focus type of equilibrium set satisfying the Shilnikov 
condition, then there are countably infinite number of periodic solutions in the 
vicinity of this homoclinic orbit. The analysis also shows that in the vicinity 
of this homoclinic orbit, complex bifurcation sequence can be expected in the 
phase portrait [|17],|l5 • 
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Another approach to explain the complex bifurcation sequences has been 
through Gavrilov - Shilnikov scenario. It has been shown that systems having 
homoclinic tangencies to periodic solutions possess quasi-random dynamics 



and MMO like behavior in the control parameter space Each of these sce- 
narios are characterized by bifurcation diagrams obtained from the stability 
analysis of the homoclinic orbit and by the corresponding scaling relations in- 
volved in the approach to the homoclinicity. Apart from these studies on the 
homoclinic bifurcations of the continuous time systems, attempts have also 
been made to study the MMOs through discrete maps [p^-[2ll . 



In order to understand the numerically obtained Poincare maps from such 
model systems and those obtained from experiments, attempts have been made 
to derive map structure starting from local analysis [|l8| , [22| , |23l . Such attempts 
have been reasonably successful in the sense that the features predicted by 
the derived maps agree with the features of the numerically obtained Poincare 
maps. However, due to the very nature of global bifurcations, there is no easy 
way of classifying the entire complex bifurcation sequences of these systems. 
Hence numerical evidence plays a crucial role in classifying the dynamics of 
these systems. 



Yet another approach to the understanding of these complex bifurcation se- 
quences and the MMOs is based on the analysis of the structure of the slow 
manifold P|, |I^ , |IT| , |26| -2B| of the system of equations. A standard slow manifold 
structure that has been used in the study of MMOs is the S— shaped struc- 
ture] 29,|30| wherein the upper and lower pleats are attractive and the middle 
branch is repulsive. By appropriately locating the fixed point on the upper or 
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lower pleat, the origin of MMOs has been explained||26||. In this construction of 
slow manifold, Shilnikov's criterion is satisfied if the direction of approach of 
the fast variable is transverse to the slow manifold containing the fixed point. 
However, in many situations |^,fn],|n|, including the system under study, the 
approach to the saddle fixed point by the fast variable is not transversal. De- 
pending on the nature of the approach of the fast variable to the fixed point, 
a classification for the MMOs has been suggested by Koper et. al. [|TJJ as 



type-I and type-II corresponding to the tangential approach and transversal 
approach (to the slow manifold) respectively. It has also been suggested that 
the occurrence of MMOs and incomplete approach to homoclinicity is related 
to the presence of Hopf bifurcation close to the fold of the S-shaped slow man- 
ifold structure [ffOl. 



In this paper, we analyze the complex dynamical behavior of a model which 
has been introduced in the context of a type of plastic instability called the 
Portevin - Le Chatelier effect. The model exhibits a rich variety of dynamics 
such as period bubbling, doubling and the complex bifurcation sequences. Two 
distinct features of the model are the atypical nature of the relaxation oscil- 
lation and the MMOs. The latter exhibits partial features of Shilnikov, and 
Gavrilov-Shilnikov scenario and shows a incomplete approach to the homo- 
clinic point. Our effort is focussed in understanding this issue in the context 
of our model. To begin with we study the nature and the origin of the re- 
laxation oscillations by analyzing the geometry of the slow manifold which 
controls the relaxation oscillations. We show that the underlying cause of the 
relaxation oscillations is due to the atypical bent geometry of the slow mani- 
fold. This feature forms the basis of further analysis of the nature of the MMO 
sequences, and the incomplete approach to homoclinic bifurcation. The paper 
is organized as follows. In section II, for the sake of completeness, we start with 
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a brief introduction to the phenomenon followed by a description of the model. 
Section III contains a detailed analysis of the complex bifurcation sequences 
exhibited by the model in the plane of two physically interesting parameters. 
Section IV contains a discussion on the origin of relaxation oscillations using 
the geometry of the slow manifold. Using this, we explain the origin of the 
MMOs and derive a scaling relation involving the maximum number of pe- 
riodic orbits allowed for given value of the control parameter. The analysis 
helps us to understand the cause of the incomplete approach to homoclinicity. 
We conclude the paper with discussion and conclusions in section V. 



2 A Dynamical Model for Jerky Flow 



Since the model is rooted in the area of plastic instability, for the sake of com- 
pleteness, we start with a brief introduction to the phenomenon. The Portevin- 
Le Chatelier (PLC) effect is a plastic instability manifesting when specimens 
of metallic alloys are deformed under tensile deformation. Under normal con- 
ditions, the stress-strain curve is smooth. However, repeated yield drops occur 
when the material parameters are in the regime of instability. Each of the load 



drops is related to the formation and propagation of dislocation bands |31| , |32 



The PLC effect (or the jerky flow) is seen in several metallic alloys such as 
commercial aluminium, brass, alloys of aluminium and magnesium [^TJ. The 
phenomenon is observed only in a window of strain rates and temperature. It 
is generally agreed that the microscopic cause of the instability is due to the 
interaction of dislocations with mobile point defects. This leads to the negative 
strain rate characteristic of the yield stress. The basic idea was formulated by 



Cottrell [33] few decades ago. However, this model and its extensions do not 
deal with the time dependent nature intrinsic to the phenomena. 
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The first dynamical description was attempted Ananthakrishna and cowork- 
ers several years ago |3^J5^]. The basic idea of the model is that most of the 
generic features of the PLC effect stem from nonlinear interactions between 
defect populations. The model in its original form does introduce spatial de- 
pendence of specific nature. However, further analysis of the model ignores 
the details of spatial inhomogeneous structure. The model consists of three 
types of dislocations and some transformations between them. The model has 
proved to be very successful in that it could explain most of the experimen- 
tally observed features such as the existence of bounds on the strain rate for 
the PLC effect to occur, the negative strain rate sensitivity, etc., |35|j36| , |39|] . 
Several aspects of the model has been investigated [H0H4!|. 



A few comments may be in order regarding the spatial aspect of the PLC 
effect. The nature of spatial terms that should be introduced in the descrip- 
tion of the PLC effect has been a controversial topic ^3f. However, there is 
some consensus that the double cross slip mechanism plays an important role 
in the spatial aspect. In the above model, justification for ignoring the spa- 
tial inhomogeneous structure and considering only the temporal aspects of the 
phenomenon is that the variables (dislocation densities) correspond to the col- 
lective degrees of freedom of the spatially extended systems. But, if one were 
to be interested in spatial aspects directly, we refer the reader to an improved 
version of the model[44j]. Further work to improve the model by including the 
nonlocal effects of immobile density along with the cross slip term is under 
active investigation. 



One important prediction of the model is that the phenomena should be 
chaotic in a certain regime of applied strain rate|J7|,[^,fIIJ. This prediction 
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has been verified by analyzing stress signals obtained from single and poly- 
crystalline samples Further, the number of degrees of freedom re- 
quired for a dynamical description of the phenomenon estimated from the 
analysis was found to be four or five consistent with that envisaged in the 
model. Moreover, since the physical system is spatially extended, this reduc- 
tion to few degrees of freedom does suggest that these modes correspond to 
collective degrees of freedom of the participating defects. From this point of 
view, dealing with the temporal aspect appears to be justified. Therefore, it is 
natural to investigate the chaotic behavior of the model in its own right. Some 
preliminary results on the chaotic aspects of the model have been published 
earlier 0,|9|,|5|. 



The model consists of mobile dislocations and immobile dislocations and an- 
other type which mimics the Cottrell's type, which are dislocations with clouds 
of solute atoms |3^| . Let the corresponding densities be N m , N im and iVj, re- 
spectively. The rate equations for the densities of dislocations are: 

N m = 6V m N m - (3N 2 m - (3N m N im + 1 N im 

-a m N m (1) 

N im = (3N 2 m - (3N i mN m - jN im + ai N h (2) 

N i = a m N m -a i N i . (3) 

The overdot, here, refers to the time derivative. The first term in Eq. (1) is 
the rate of production of dislocations due to cross glide with a rate constant 9. 
V m is the velocity of the mobile dislocations which in general depends on some 
power of the applied stress a a . The second and third term refer to annihilation 
or immobilization processes. The fourth term represents the remobilization of 
the immobile dislocations due to stress or thermal activation (see the loss 
term 7iVj m in Eq. 2). The last term represents the immobilization of mobile 
dislocations either due to solute atoms or due to other pinning centers. a m 
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refers to the concentration of the solute atoms which participate in slowing 
down the mobile dislocations. Once a mobile dislocation starts acquiring solute 
atoms we regard it as a new type of dislocation, namely the Cottrell's type iVj. 
This process is represented as a gain term in Eq. (3). As they acquire more 
and more solute atoms they will slow down and eventually stop the dislocation 
entirely. At this point, they are considered to have transformed to N im . This 
process has been represented by the loss term in Eq. (3) and a gain term in 
Eq. (2). 

These equations should be dynamically coupled to the machine equations de- 
scribing the rate of change of the stress developed in the sample. This is given 
by 

° a = K(e' a -B N m V m ), (4) 

where k is the effective modulus of the system, e a is the applied strain rate, 
V rn is the velocity of mobile dislocations and B is the Burgers vector. These 
equations can be cast into a dimensionless form by using the scaled variables 



x = N m (Pj,y = N im (J^, 

z = nJ ,t = 6V t, and <f> = (^) . 

Using the power law dependence V m = V (^) m , Eqs. (1-3) and (4) can be 
rewritten as 



x = (j) m x — ax — b x 2 — xy + y, (5) 

y = b (b x 2 -xy-y + az), (6) 

z — c(x — z), (7) 

= d( e -0"V), (8) 

Here a = (a m /9V ),b = (>y/8V ),c = (ai/9V ), k = (9f3a d/>yB ) and 
e = (e a (3 / BqVq^) . For these set of equations there is only one steady state 
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which is stable. There is a range of the parameters a, bo, c, d, m and e for which 
the linearized equations are unstable. In this range x, y, z and are oscillatory. 



Among these physically relevant parameters, we report here the behavior of 
the model as a function of most important parameters namely the applied 
strain rate e and the velocity exponent m. We use e as the primary control 
parameter for the analysis. The values of other parameters are kept fixed 
at a = 0.7, b = 0.002, c = 0.008, d = 0.0001 and k = 1.0. The present 
choice of parameters does not necessarily correspond to a realistic experimental 
situation, although there is a range of allowed values. As can be verified these 
equations exhibit a strong volume contraction in the four dimensional phase 
space. We note that there are widely differing time scales corresponding to 
a, bo, c and d (in the decreasing order) in the dynamics of the model. For this 
reason, the equations are stiff and numerical integration routines were designed 
specifically to solve this set of equations. We have used a variable order Taylor 
series expansion method as the basic integration technique with coefficients 
being determined using a recursive algorithm. Most of the bifurcation analyses 
were performed using these indigenous routines. AUTO software package |50[ 
was used exclusively for two parameter continuation of bifurcation points. 
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3 Summary of bifurcation exhibited the model 

In an attempt to understand the complex bifurcation sequences exhibited by 
the model, we start with an outline of the gross features of the phase diagram 
in the (m, e) plane shown in Fig. 1. In our discussion, we consider m to be the 
unfolding parameter. For values of m > ma ~ 6.8, the equilibrium fixed point 
of the system of equations is stable. At m — m d , we have a degenerate Hopf 
bifurcation as a function of e. For values less than m^, we have a back-to-back 
Hopf bifurcation. The periodic orbit connecting these two Hopf bifurcations is 
referred to as the Principal Periodic Orbit (PPO). The dynamics of the system 
is essentially bounded by these two Hopf bifurcations. In Fig. 1, the broken 
line represents the Hopf bifurcation and the dotted lines correspond to the 
first three successive period doubling bifurcations leading to period 2, 4 and 
8 orbits. The region in between the first period two and the Hopf bifurcation 
line exhibits monoperiodic relaxation like oscillations. 

The PPO for most of the parameter plane (m, e) is born through a subcritical 
Hopf bifurcation leading to relaxation like oscillations. However, the narrow 
region between the Hopf bifurcation line (corresponding to large values of e) 
and the period two regime is characterized by small amplitude nearly symmet- 
ric coplanar limit cycles. The complex bifurcation sequences, characterized by 
alternate periodic-chaotic sequences seen in the parameter space, are roughly 
indicated by the hatched region. Since the part of the hatched region extends 
beyond the outermost period doubling line (large values of e and small values 
of m), both the MMOs and the small amplitude monoperiodic limit cycle so- 
lutions coexist in this region. (See also Fig. 3.) A codimension two bifurcation 
point in the form of a cusp (shown as a filled diamond) at (e c , m c ) formed by 
merging of the locus of two saddle node periodic orbits (represented by bold 
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lines) of the PPO is also displayed in Fig. 1. Apart from these bifurcations, 
we failed to detect any other bifurcation or equilibrium set in the phase space. 
We will deal with each of these regions in detail below. 



Bifurcation diagrams have been obtained by plotting the maxima of any one of 
the variables x,y, z or as a function of the control parameters (e, m). We have 
mostly shown the bifurcation diagrams in the variable x. This choice enables 
a good visual representation of the bifurcation diagram since the maxima of 
the x variable is quite large compared to other variables. Based on the nature 
of the bifurcation sequences, the parameter space can be broadly grouped into 
two regions, viz. m > 2.0 and m < 2.0, and we will discuss the changes as a 
function of e fixing m at a particular value. 



3.1 Region m> 2.0 



We briefly summarize the results for this region. For values of m > 2.16, 
the bifurcation diagrams are characterized by an incomplete period doubling 
cascades followed by reverse period doubling bifurcations, displayed as nested 
bubbles of periodic states (see Fig. 2). As m decreases, the number of periodic 
bubbles nested in the structure increases as 2 n , with n — > oo, culminating in 
chaos. Just below m = 2.16, the disjoint chaotic bubbles collide with each 
other forming an extended attractor which has been referred to as 'bubble 
bursting' in the literature [^TJ. Similar features have been observed with m as 
the control parameter keeping e fixed at an appropriate value. The rates of 
the period doubling (PD) bifurcations as well as the reverse period doubling 
bifurcations with respect to e and m fall close to the value of the Feigenbaum's 
constant for quadratic unimodal maps, namely b~p = 4.66. When the value of 
m reaches a critical value m = 2.11, a period three cycle is born through a 
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saddle node bifurcation and has the largest width (in e) for m ~ 2.0 in this 
regime. 

3.2 Region m < 2.0 

For m < 2.0, the system exhibits qualitatively different behavior compared 
torn > 2.0, in the sense that higher order MMOs emerge gradually as m is 
decreased. In this region, the system exhibits complex bifurcation sequences 
or the alternate periodic-chaotic sequences which are characteristic to this sys- 
tem. The stable periodic orbits in the bifurcation sequence typically exhibit 
MMO nature and they are labelled by L s , where L is the number of large 
amplitude loops and s is the number of small amplitude loops of the periodic 
orbit. These MMOs are heralded by the creation of the period three (L s : l 2 ) 
region after the PD cascade to chaos in the bifurcation plot. To illustrate the 
nature of bifurcations in this region, we fix the unfolding parameter at m = 1.8 
and discuss the bifurcation with respect to e. 

Figure 3 shows the bifurcation sequence with alternating chaotic and periodic 
states along with the higher order periodic isolas (isolated bifurcation curves). 
The unstable periodic orbits are shown by dashed lines. In the case of isolas, 
we have shown only the largest amplitude isola for any given periodicity. The 
first instability of the PPO through a PD bifurcation opens up a period dou- 
bled orbit having a large parameter width in e. This feature persists for the 
entire m c < m < 2.0 regime. As in the case of m — 2.0, a period three isola 
is born through a saddle node (SN) bifurcation from the chaotic attractor. 
In Fig. 3, stable periodic orbits are shown to be bounded between PD bifur- 
cation (filled circle) and the SN bifurcation (filled triangle). The sequential 
way SN and PD bifurcations arrange themselves to form an isola can be easily 
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understood by considering the behavior of the Floquet eigenvalue of the pe- 
riodic orbit. The period three orbit is born in a SN bifurcation accompanied 
by the disappearance of the first chaotic window attractor. At this value of 
e, Floquet eigenvalue is at +1 creating a pair of stable and unstable period 
three orbits. As e is increased, the eigenvalue of the unstable periodic orbit 
increases beyond +1 while the eigenvalue corresponding to the stable orbit 
keeps decreasing and crosses —1 resulting in a PD bifurcation. Due to the 
isola structure, further increase in e makes this eigenvalue cross back —1 thus 
restabilizing the period three orbit. For any further increase in e, the Floquet 
eigen values of the stable and unstable orbit merge again at +1 and vanish 
in another SN bifurcation to complete the isola structure. The next higher 
order isola is also created through a SN resulting in the disappearance of the 
chaotic attractor born from the destabilisation of stable period three orbit. 
Higher order periodic orbits (isolas) are formed in a similar way. Note that 
these isolas are independent of the PPO. We refer to this sequence of periodic 
orbits of this form as the principal period adding sequence (PPAS) or the 
principal Farey sequence. In the MMO notation, the PPAS can be written as 
U : l n where n = 2,3,4,.... As the periodicity increases, the width of the 
isolas in e decreases. Since the isolas are independent of the PPO, any change 
in the stability of the PPO has no effect on the nature of the sequence. This is 
evident from the bifurcation diagram where higher periodic orbits (isolas) are 
seen even after stability of the PPO is reestablished. The restabilisation of the 
PPO is through a reverse period doubling cascade from chaos. This chaotic 
segment formed from the reverse period doubling of the PPO expands in an 
interior crisis due to its collision with an unstable periodic orbit of the next 
higher period isola (period three in the case of m — 1.8). This crisis point, 
shown by an arrow, marks the lower boundary of the multistability region in 
the bifurcation diagram. 
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The inset of Fig. 1 shows the expanded region of interest of the phase dia- 
gram in the (m, e) plane. As can be seen, the locus of the SN bifurcations 
corresponding to the MMOs are distinct and higher period SN bifurcation 
curves cross the lower ones resulting in the period n isola extending beyond 
the period (n — 1) isola. Moreover, the region where SN bifurcation curves 
overlap with the region of the period doubling curves, multistability regions 
in the parameter space (m, e) are created. This is clearly seen in Fig. 3 where 
bifurcation diagram for m — 1.8 is shown. (See also Fig. 4.) Typically, the 
same mechanism as described for the case of m — 1.8 operates for the period 
adding sequences in the region 2.0 > m > m c , where m c ~ 1.1 is the value of 
m at the cusp point. As m decreases from m = 2.0, higher number of stable 
periodic windows are accommodated with concomitant decrease in the width 
of the chaotic regimes separating the periodic windows. The arithmetically 
increasing periods of the orbits going from left to right form an incomplete 
period adding sequence with decreasing widths for higher order periodic win- 
dows. These features are shown in the bifurcation diagrams for m = 1.8, and 
m — 1.2 in Fig. 3 and Fig. 4 respectively. 

Below m c , the bifurcation diagram is even more interesting and rich. Here, we 
outline the features related to the Farey states. For the case m = 1.0, only 
three principal Farey states denoted by L s , s = 1,2 and 3 survive, as shown 
in Fig. 5. The well developed sub- Farey sequences are also shown in the inset 
of Fig. 5. The sub-Farey states created go from right to left in contrast to the 
principal Farey states (see Fig. 5). All these sub-Farey sequences culminate in 
a SN bifurcation. While in the first bifurcation of the PPO (SN bifurcation), 
the transition is from 1° — > oo 1 , the mid region of parameter accommodates 
both the large amplitude and small amplitude solutions with nearly equal 
weights. Towards the end region of e, we find no fine structure typical of the 
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first two principal chaotic windows. 
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4 Mechanism of relaxation oscillations 



One characteristic feature of the dynamics of the system is its strong relaxation 
nature. This is seen even in the case of the mono-periodic solutions emerging 
from the Hopf bifurcation for small values of e. This feature, of course persists 
for other regions of the (m, e) plane where the MMO type of oscillations are 
also seen. These two aspects are interrelated and are a result of the structure of 
the slow manifold as we will show. Since our system does not follow the known 
homoclinic scenarios, we look for a new mechanism for the MMOs based on 
the mechanism proposed for the relaxation oscillations. 



4-1 Relaxation oscillations 



Relaxation oscillations are highly nonlinear oscillations with large amplitude 
excursions of the fast variable. These oscillations arise as a consequence of the 
existence of a fast time scales compared to the time scales of other variables 
in the dynamics of the system. The relaxation oscillations have been an in- 
tense area of research in the context of biological rhythms [|52]. The relaxation 
oscillations that manifest in the model under study is a type of relaxation os- 
cillation wherein the fast variable takes on large values for a short time after 
which it assumes small values of the same order of magnitude as that of the 
slow variables. The time spent by the fast variable in the part of phase space 
where the amplitude of the fast variable is small is a substantial portion of 
the period of the orbit. It is this type of relaxation oscillation that is domi- 
nantly seen in our system, even though other types of relaxation oscillations 



are also seen [36] for certain other regimes of the (m, e) plane. We shall refer 
to this type of relaxation oscillation as pulsed type relaxation (PTR) oscilla- 
tion. A typical plot of x(t) is shown in inset of Fig. 6 for e = 200.0 and m = 1.2. 
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In the case when two disparate time scales are present in the dynamics, using 
multiple scale perturbative analysis, Baer and Ernaux have shown that Hopf 
bifurcation can lead to relaxation type of oscillations |53|. They have shown 
that nearly sinusoidal solutions born out of the Hopf bifurcation change over 
to relaxation like oscillations in a small region of the value of the control pa- 
rameter. In such the crossover to relaxation oscillations is confined to 
the slow manifold around the fixed point. As we will see, the nature of the re- 
laxation oscillation in our model is very different from that discussed by Baer 
and Ernaux. Here, it suffices to say that the PTR is a result of the evolution 
wherein the trajectories visiting the slow manifold region around the fixed 
point are pushed out to another part of the slow manifold away from the fixed 
point where trajectories spend considerable fraction of its period. 



To understand the nature of the relaxation oscillations, we first study the 
structure of the slow manifold (S) and the behavior of the trajectories visiting 
different regions of S. Consider the slow manifold given by 

x = g(x, y, 0) = —b x 2 + xS + y = (9) 

with 5 = (f) m — y — a. Here, the slow variables y and <fi (and therefore 5) are 
regarded as parameters. Further, as we will see below, it is simpler to deal 
with the structure of the slow manifold in terms of the S instead of both y 
and (p. Then, the physically allowed solution of the above equation is 

5 + VP + Ab y 
X = 2b (10) 

where 6 can take on both positive and negative values. Noting that &o is small 
and therefore 5 2 ^> 4boy, two distinct cases arise corresponding to 5 > and 
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5 < for which x ~ 5/b and a; ~ — y/5 respectively. Further, since the slow 
variable and y take on values of the order of unity, the range of 5 = 5(y, 0) 
is of the same order as that of <fi and y (as is evident from Figs. 6 and 7). 
Thus, we see that x ~ —y/5 is small and x ~ 5/b is large. For values around 
5 = and positive, we get x ~ (y/bo) 1 ^ 2 - 



Let Si denote the region of slow manifold values of x corresponding to 5 > 
and S 2 the region of slowmanifold values of x for 5 < 0. The bent-slow mani- 
fold structure along with the two portions of the slow manifold are shown by 
a bold lines in the (x, 5) plane in Fig. 6. A local stability analysis for points 
on S\ and S 2 shows that dg/dx = 5 — 2b Q x is negative. This implies that 
the rate of growth of x is damped and hence these regions, Si and S 2 will be 
referred to as attracting. In Fig. 6, we have shown a trajectory corresponding 
to a mono-periodic relaxation oscillation (m = 1.2 and e = 200.0) by a thin 
line. As can be seen, the trajectory spends most of the time on Si and S 2 . 
For points below the line 2b x = 5 (5 > 0), dg/dx > implying a positive 
rate of growth of the x variable and hence we call this region as repulsive or 
'unstable' (shaded region of Fig. 6). We stress here that this region is not a 
part of the slow manifold. Even then, the trajectory starting on S 2 does con- 
tinue in the direction of increasing 5 beyond 5 = 0. Once the trajectory is in 
this region, it moves up rapidly in the x direction (due to the 'unstable' na- 
ture) until it reaches x = 5/2b Q line, thereafter, the trajectory quickly settles 
down on to the Si part of the slow manifold due to the fact dg/dx becomes 
negative. As the trajectory descends on Si approaching S 2 , we see that the 
trajectory deviates away from Si. This happens when the value of x is such 
that 2box < 5, i.e., dg/dx > 0. Thus, points on Si satisfying this condition are 
locally unstable. (For points in this neighborhood 5 ~ 0.2, a; ~ 50.) Thus, the 
trajectory makes a jump from Si to S 2 in a short time. This roughly explains 
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the origin of the relaxation oscillation in terms of the reduced variables 5 and x. 

The actual dynamics is in a higher dimensional space and a proper under- 
standing will involve the analysis of the movement of the trajectory in the 
appropriate space. Moreover, quite unlike the standard S— shaped manifold 
with upper and lower attracting pleats with the repulsive (unstable) branch, in 
our model, both branches of the bent-slow manifold are connected, and there 
is no repulsive branch of the slow manifold. Thus, the mechanism of jumping 
of the orbit from S 2 to S\ is not clear. In order to understand this, consider 
a 3-d plot of the trajectory shown in Fig. 7. Retaining the same notation for 
the 3-d regions of the slow manifold as that used for the x — 5 plane, regions 
Si and S 2 are shown in Fig. 7. As can be seen, the region S 2 corresponding 
to small values of x lies more or less on the y — plane and the region S± 
corresponding to large values of x is nearly normal to the y — 4> plane due to 
the large 6g 1 factor. (Note that the scales of y and are the smallest for the 
system. ) Regions Si and S 2 are demarcated by the 'fold curve' which lies in 
the y — <p plane and is given by 5 = <p m — y — a = 0. As in the case of x — 5 
plane, in 3-d space also, the rapidly growing nature of the trajectory seen in 
the approximate region below the surface of 2box = (jf 1 — y — a and lying to 
right of the 'fold curve' is due to dg/dx > 0. 

The principal features of the relaxation oscillations that we need to explain 
are: a) very slow time scale for evolution on S 2 , b) fast transition from S 2 to 
Si and c) evolution on Si. As mentioned in the introduction, these are related 
to the slow - fast time scales. In order to understand this, we shall analyze 
Eqs. (6) and (8) by recasting them in terms of 5. In the whole analysis it 
would be helpful to keep in mind the range of values of x, y, z and (shown 
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in Figs. 6 and 7), in particular, their values as the trajectory enters and leaves 
Si. Consider rewriting Eq. (6) valid on the slow manifold S in terms of 5: 



y = b (x5 — xy + az). 



(11) 



The idea is to study this equation along with Eq. (8) in specific regions of the 
phase space to understand the general features of the flow, viz., on S2, just 
outside S 2 , and on Si. The presence of the z variable in Eq. (11) poses some 
problems. However, it is possible to get a rough estimate of the magnitude of 
z and the relative changes in the values of z which is all that will be needed for 
our further discussions. To see this, consider Eq. (7) from which we see that z 
follows x. Further, Once the trajectory moves out of S2, x changes rapidly and 
therefore the value of z increases (in a relatively short interval of time), reach- 
ing its maximum value, z max , just before the trajectory returns to S 2 . When 
the trajectory is on S 2 , since x ~ y/\S\, from Eq. (7) we see that the value of z 
is slowly decreasing (with a time constant c _1 ) starting from z max , reaching its 
minimum value, say z m i n , around the time when the trajectory leaves S 2 . In 
other words, the magnitude of z is maximum when the trajectory enters S 2 and 
minimum when it leaves S 2 . Further, we note that x = —y/S <C x = z ~ e/2 
and z oscillates around its equilibrium value zq. Thus, the values of z max and 
z m in are larger than the range of allowed values of y. (Note that this is also 
consistent with the fact that the time scale of z is larger than that for y and (f>.) 

Consider the behavior of Eq. (11) on S 2 . Using the values of x ~ y/\S\, we get 



By noting that on S 2 , z decreases from z max to z m i n , we see that there is 
a range of small values of y for which y > and for relatively larger values 




+ az 



(12) 
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of y, y < 0. Thus, y clearly has a turning point on S2 beyond which y decreases. 

Next, consider the changes in 0. Using the value of x — y/\5\ on S 2 in Eq. 
(8), we find that e is much larger than <f) m y/\5\, since these variables are of 
the order of unity. Thus, increases linearly, at a rate close to de « 1. Con- 
sidering the fact that x ~ for the entire interval the trajectory is on S 2 , the 
time scale of evolution of the trajectory is entirely controlled by the two slow 
time scales of y and <fi. This roughly explains the behavior of the trajectory 
on S 2 - 

Now consider the behavior of Eq. (11) in a small region just outside 5 = 0. 
Using x ~ (y/bo) 1 ^ 2 valid for 5 > 0, we get, 



We are interested in investigating the behavior of Eq. (13), for z = z max , ap- 
propriate as the trajectory approaches S 2 and z = z m i n , appropriate as the 
trajectory leaves S 2 . Consider the first choice (z = z max ) corresponding to the 
trajectory as it approaching S 2 from outside (5 > 0). Then, using the value of 
b , an order of magnitude calculation shows that there is a range of small val- 
ues of y for which y > 0. This implies that y grows for small y, meaning that 
the trajectory moves towards S2- Now consider using z = z min corresponding 
to the situation when trajectory has left S2. Similar estimation shows that 
there is a range of (relatively larger) values of y for which y < 0. This implies 
that y decreases for relatively large values of y, meaning that the trajectory 
is moving away from S2. ( Note that for this case, there may or may not be 
a range of y for which y > 0.) Thus, in both cases the directions of growth of 
y for small and large y just outside S2 are consistent with the behavior of y 




(13) 
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just inside S 2 . (See Fig. 7.) 



Now, consider Eq. (8) with x ~ (y/bo) valid for the region 5 positive but 
small. Then, 

(14) 

Keeping in mind the order of magnitude of bo, and the fact that y and are 
of the order of unity, the magnitude of <p m (y/bo) 1 ^ 2 is seen to be larger than 
its value on S 2 - Note that a quick order of magnitude calculation shows that 
there are values of y and such that m (y/bo) 1 ^ 2 is of the order of e which 
implies that is about to decrease and therefore is near its maximum. More- 
over, if anything, <J) m x in Eq. (8) increases as the trajectory tends to moves 
out of S 2 , since x ~ xS just outside S 2 . This implies will eventually decrease. 

Combining the results on y and for regions just outside and inside the 'fold', 
we see that the trajectory enters S2 in the region corresponding to small values 
of y and 0, and makes an exit for relatively larger values of and y (compared 
to their values as the trajectory enters S 2 ). Finally, we can see that just to 
the right of 5 = line, x ~ x5, with 5 very small, which suggests that the 
time constant is small. Thus, the growth of x is slow in the neighborhood of 
5 = 0, and is tangential to the S2 plane even in the 'unstable' region. However, 
once the trajectory moves away from 5 = 0, the growth of the trajectory is 
controlled by dg/dx and hence the time scale of growth of x is of the order of 
S^ 1 which is of the order of unity. This essentially explains why the trajectory 
tends to move into the 'unstable' region and grows rapidly. 

Once in the 'unstable' region, the value of x continues to grow in this region 



(f) = d 



e — 



bo 



1/2 
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of the phase space as can be seen from Eq. (8) until the value of x is such that 
(f) m x = e is satisfied. Beyond this value of 0, <fi is negative. Thus, the trajectory 
leaving S2 eventually falls onto the Si part of the slow manifold. We can again 
evaluate y and <j) just as the trajectory reaches S±. Using x ~ 5/b in Eq. (11), 
we find 

(15) 

The sign of y is determined by the factor (5 — y) at the point where the 
trajectory reaches Si. To see the relative magnitudes of 5 and y, consider 
obtaining an equation for 5 starting from g(x,y,5) = 0. Differentiating this 
and using x = on for the slow manifold, we get 

S = -—[x(8-y) + az]. (16) 



y = b 



— (5-y) + az 



Using x ~ 5 /bo on Si, we see that y is a fast variable compared to 5. Thus, in 
this interval of time, we could take y — 0, i.e., 



bnaz 



(17) 



Since all these variables 5, y and z are positive on Si, we see that y > 5. (Note 
the factor boaz/5 is small.) Using this in Eq. (15) we see that y decreases. 
Now, consider the equation for <f>. Using x ~ S/b on Si, we get 

<f) m 5~ 



cf) = d 



bo 



(18) 



Noting the value of bo, we see that <fi will be negative when the trajectory 
reaches Si. The time scale of evolution of y in Eq. (15) is of the order of 
unity while that of is ~ d/b . These time scales are relatively fast. (These 
statements are true only as the trajectory hits Si.) Moreover, since x is a 
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fast variable, the changes in x component dominates the descent of the tra- 
jectory. Finally, as the trajectory approaches S 2 , dg/dx becomes positive and 
the trajectory jumps from Si to S 2 . Combining these results, we see that the 
trajectory moves towards the region of smaller values of y and (f> entering S 2 
in a region of small values of y and 0. 

In summary, the sequential way the orbit visits various parts of the phase 
space is as follows. The trajectory enters £2 part of the slow manifold in re- 
gions of small y and <fi making an exit along S 2 for relatively large and 
y. Thereafter, the trajectory moves through the 'unstable' part of the phase 
space before falling onto the Si and quickly descends on Si. This completes 
the cyclic movement of the trajectory and explains the geometrical feature 
of the trajectory shuttling between these two parts of the manifold and the 
associated time scales. 

Now, the question that remains to be answered is — do the trajectories always 
visit both Si and S 2 or is there a possibility that the trajectory remains con- 
fined to Si ? It is clear that if the former is true, relaxation oscillations with 
large amplitude will occur and if the latter is true, the oscillations are likely to 
be of small amplitude. Here, we recall that the coordinates of the saddle focus 
fixed point are xq = z ~ e/2 which is much larger than the values of x on S 2 
(~ yj Thus, the fixed point located on the Si will be close to the 'fold' 
at the first Hopf bifurcation, e = e/, since the latter occurs at small values 
of e (ey ~ 5). Due to the unstable nature of the fixed point, the trajectories 
spiralling out are forced onto the S 2 part of the manifold resulting in relax- 
ation oscillation. This point has been illustrated by considering the example 
of a period eleven orbit for m — 1.2 and e = 267.0 shown in Fig. 8. As is clear 
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from this diagram, the small amplitude oscillations are located on the Si. As 
the small amplitude oscillations grow, the relaxation nature does not manifest 
until the orbit crosses over to S2. To the best of the authors knowledge, the 
mechanism suggested here for pulsed type relaxation oscillations is new. 



The above feature of the trajectories continuing in the same direction of the 
slow manifold (S2) well into the 'unstable' part of the phase space is some- 
what similar to canard solutions where the trajectories tend to follow the slow 
manifold well into the repulsive part of the slow manifold before jumping to a 



attracting branch [p|p4|. The differences, however, are clear. While in canard 
solutions, the trajectory tends to move along the repulsive part of the slow 
manifold before jumping to the attracting branch of the slow manifold, in our 
case, the trajectory leaves the slow manifold and moves into the 'unstable' 
part of the phase space which is not a part of the slow manifold. 

4-2 Mixed Mode Oscillations 

We now consider the origin of the MMO sequences in our model. Global bi- 
furcation scenarios are known to be relevant to the MMOs and in the intro- 
duction, we briefly mentioned two of the possible global bifurcation scenarios 
which display MMO like sequences. These scenarios are based on the homo- 
clinic contacts of an equilibrium set like the saddle focus fixed point and saddle 
periodic orbit for the Shilnikov and Gavrilov-Shilnikov scenarios respectively. 
Each of these scenarios are characterized by the bifurcation diagrams obtained 
from the stability analysis of the homoclinic orbits and by the corresponding 



scaling relations involved in the approach to homoclinicity [17,23 . 



First, we consider the similarities of the behavior of our model with the char- 
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acteristic features of the Shilnikov scenario. In three dimension, the Shilnikov 
criterion is stated in terms of the two possible combinations of the dimen- 
sions for the invariant manifolds of the saddle-focus; the unstable manifold is 
two dimensional and stable manifold is one dimensional and vice versa. For 
these two cases, the presence of an homoclinic orbit is given by the condition 
|p/A| < 1, where the eigen values of the fixed point are given by p ± iu>, —A, 
where p > 0, A > for the first case and p < 0, A < for the latter. The analy- 
sis of the Shilnikov scenario shows that in the neighborhood of the homoclinic 
point, the parameter space is organized such that the period of the principal 
periodic orbit tends to infinity as the parameter approaches the value corre- 
sponding to the homoclinic point. In our case, the system is four dimensional, 
with the unstable manifold of the fixed point characterized by a pair of com- 
plex eigenvalues p ± iu (p > 0) and the stable manifold by two eigenvalues 
Ai < and A2 < 0. Here, Ai stays close to zero and A2 is substantially neg- 
ative. Thus, in our case, the criterion \p/\\ < 1 refers to \p/\\\ < 1. We find 
that this condition is satisfied only in a small region just prior to the disap- 
pearance of the PPO in a Hopf bifurcation. A typical plot of the eigenvalues 
for m — 1.2 is shown in Fig. 9, where we have also shown the phase uj. Even 
though \p/\\ < 1 is not satisfied over large portion of e and m, we do see that 
the period of the periodic orbits tend to increase as e is increased (m < 2.0) 
which is typical of the Shilnikov scenario. A plot showing the period (of the 
superstable orbits) verses the deviation from estimated homoclinic point(e*) 
is displayed in Fig. 10 for m = 1.4. (We have plotted points from period three 
onwards.) Here, we have taken the value of e* to be the value of e for the onset 
of the last observed periodic orbit with period 12 (e = 247.63). It must be 
stated that we do not face any difficulty in locating any of the periodic orbits 
upto the period 12. However, we are unable to detect the next period which 
we interpret as an incomplete approach to the true homoclinic point. Here, it 
must be mentioned that incomplete approach to homoclinicity is quite com- 
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mon [^^^,|6|j8| JT0|JTT|JT3|j27f] . We stress that even in the region where Shilnikov 
criterion is obeyed, we do not observe homoclinic orbit. 



One other feature which is usually seen in the Shilnikov scenario is that the 
reinjection of the trajectory in to the neighborhood of the fixed point is along 
the direction of the fast variable after which the trajectories tend to stay 
around the saddle focus fixed point. In our model, since x is the fastest vari- 
able, it also acts as a reinjection direction. However, a closer examination 
shows that the spiraling in of the orbit towards the fixed point is along the z 
direction which is the next fastest variable. This is evident in Fig. 11 where 
a typical trajectory is shown. Even more dominant feature of the Shilnikov 
scenario is that the successive bifurcations should be connected by the PPO. 
This, however, is not true in our case as we have seen earlier, since the isolas 
which form the period adding sequence are distinct from the PPO. This fea- 
ture is clear from Fig. 3 for m — 1.8. In fact, this feature of the isolas being 
distinct from the PPO is more like that of the Gavrilov- Shilnikov scenario 
which requires the presence of an unstable periodic orbit which we failed to 
detect in the entire parameter region wherein nontrivial dynamics is present. 
This rules out the possible presence of any homoclinic bifurcation due to the 
saddle periodic orbit. Thus, we see that our model has partial features of both 
these scenarios. 

The above discussion suggests that the origin of the MMOs in our model is 
likely to be different from the two scenarios. In order to understand the mech- 
anism of the MMOs in our model, we will use the information on the nature 
of the relaxation oscillations. We first note that the fixed point (x ,y , z , <fi ) 
is on the S\ part of the slow manifold and moves up on Si as e is increased. 
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Both x and z have a near linear dependence on e namely xq = Zq ~ | , while 
and 0o are practically constant. Since the fixed point is unstable, any orbit 
in its neighborhood will locally expand along its unstable directions. Thus, we 
expect to get insight into the mechanism of the MMOs by studying the rate 
of expansion of such orbits. In order to understand the mechanism operating 
in our model, let us consider a periodic orbit of L s : l 10 type shown in Fig. 12. 
If the orbit has reached the neighborhood of the fixed point, any orbit on Si 
should spiral out with a local dynamics determined by the linearized eigen- 
values around the fixed point. Within this approximation, the orbit expands 
at rate exp[2np/uj] per rotation around the fixed point. Assuming a linearized 
behavior for n rotations ( for fixed values of e and m), we get 

r n /n = exp[2n(p/u)(n-l)] (19) 

where r n is the distance measured from the fixed point after n rotations along 
a fixed direction in the unstable manifold of the fixed point. Here, r 1 is the 
value of r n for n — \. Since Eq. (19) is based on linearized approximation, the 
values of r n obtained from the phase plots will be in general different due to 
the influence of nonlinearities. For this reason, we will first study the region 
of validity of Eq. (19). We note that the unstable manifold is nearly in x — z 
plane and in the neighborhood of the fixed point, the major contribution to r n 
comes from x and z. Thus, it would be sufficient to consider x(t) (or z{t)) in 
place of r n and in particular, we will use the minima or maxima to analyze the 
small amplitude oscillations of the periodic orbit using Eq. (19). Consider the 
plot shown in Fig. 12. We note that the time interval between the successive 
minima or maxima of x(t) can be taken to correspond to one rotation of the 
orbit. We shall denote the the deviations of the n-th minimum from the fixed 
point value, x (e) — x™ m , by x* n . In Fig. 12, x(t) and z(t) are plotted along with 
their fixed point values (x = z = 137.82). ( We have also shown regions of 
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x(t) corresponding to regions Si and S 2 of the slow manifold.) The decreasing 
nature of the maxima values of the amplitude of x(t) for the first few cycles 
is a reflection of the fact that the orbit has not reached the neighborhood of 
the unstable manifold of fixed point lying on Si (see Fig. 11). It is only later 
(in time) that the expanding nature of the oscillations manifest, seen as the 
increase in the magnitude of the successive maxima values of x(t). Thus, the 
value of x\ read off from x(t) (in place of r\) will contain contributions arising 
from reinjection mechanism. Hence, identifying x\ as representing the value 
of the first minimum of x(t) would be an incorrect, if one wishes to use Eq. 
(19). Thus, x\ has to be estimated by extrapolating the values of x* using 
values of n where nonlinearity plays an insignificant role, i.e., n > 3 for the 
case of Fig. 12. We denote this extrapolated value by x\. In addition, as the 
amplitude grows as a function of n, one should also expect that the linear dy- 
namics breaks down. Thus, for larger n values, we should again see the effect 
of nonlinearity. In fact, this feature shows up as a decrease (though marginal) 
in the time interval between successive minima for higher n as can be verified 
from Fig. 12 (seen between the ninth and the tenth peaks). 

Now, we attempt to estimate the changes in the magnitude of the small am- 
plitude oscillations located on Si as a function of e and estimate at what 
value of e the trajectory hits the 'fold' between Si and S 2 . We denote the 
distance of the fixed point (xo(e), yo(e), zo(e), <f>o(e)) from the 'fold' given by 
D = ((x o (e)- a ;o(e / )) 2 + (i/ o (e)-yo(e / )) 2 + (2 ; o(e)-2 ; o(e / )) 2 + (0 o (e)-0o(e / )) 2 )i 
where (xo(e/), yo(ef), zo(ef), (/>o(e/)) refers to the value of the fixed point at the 
first Hopf bifurcation (e = e/). Further, we note that the fixed point is close to 
the 'fold' at the first Hopf bifurcation and noting x = z ~ | with a very weak 
dependence of y and O on e - Thus, in one dimension where we are dealing 
with x variable alone, we can take D ~ xo(e) — xo(ef). (Note that position of 
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'fold' is insensitive to e.) Using the fact that x (ef) is small, we get D ~ x (e). 
Using this, we will attempt to find the maximum value of n allowed for which 
the condition x* n > D is satisfied. For this, we need to know the dependence 
of p and wone which decides the rate of growth of the small amplitude os- 
cillations as a function of e. In the range of the MMO sequences that we are 
interested, p decreases and ui increases nearly linearly as can be seen from Fig. 
9. Thus, we take p = p — m p e and uj = u + m w e, where, m p and m w are the 
corresponding slopes. These are evaluated numerically as the best fit for the 
region of interest of e. The fit yields p = 0.10433, m p = 0.0003632, lu ks and 
= 0.0004114 for the case shown in Fig. 12. Using this in Eq. (19) gives 

4 = exp[27r(p/^)(n - 1)] = exp^ + k 2 /e](n - 1)} (20) 
x 1 

Here k\ and k 2 are functions of p , and m p . This equation can be inter- 
preted as a scaling form for the small amplitude oscillations of stable periodic 
orbits (i.e., for a fixed n) as a function of e. Now, consider the set of all sta- 
ble periodic orbits of the form L s for a given value of m. Then, for each of 
these L s orbits, the values of n ranges upto s. Since the magnitude of small 
amplitude oscillations of these periodic orbits depends e, we can plot lnx* n 
versus 1/e, where we have used x* n values corresponding to the n-th minimum 
of a periodic orbit. The plots of In function of e 1 for n = 4 to 7 are 

shown in Fig. 13. In the figure, the lowest curve (+) corresponds to n = 2 
and other successive higher curves refer to n = 3 upwards. The lowest band 
within the n = 2 curve corresponds to a periodic orbit of the form L w and 
successive bands have decreasing s values. (The gaps correspond to chaotic 
bands between successive periodic orbits.) It is clear that the plots are linear 
and the slopes of the curves corresponding to n = 4, 5, 6 and 7 show an in- 
creasing trend, increasing in multiples of 1530 (k 2 ). To illustrate the presence 
of nonlinearity, we have also shown plots of 1m* for n = 2 and 3. One can 
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easily notice that these two curves deviate from linearity considerably for large 
e. In addition, we see that the slopes of these two lines are not in multiples 
of &2- This suggests that it should be possible to collapse all the curves for 
n > 3 onto a single curve. Noting that the slopes of curves are in multiples of 
k 2 , lnx*Jn renders them parallel and the preexponential factor which satisfies 
Eq.(20), denoted by x\ can be easily determined. The value of x\ so obtained 
can be now used in 

4 = exp[27r(p/c)(n-l)] (21) 
x[ 

to estimate the maximum number of small amplitude oscillations allowed for 
any given value of e before the size of the orbit (under the linear approxi- 
mation) is large enough to hit the 'fold'. This is determined by the value of 
n = n c at which x* n = x* c > D ~ e/2. 

In the inset of Fig. 13, we have shown that the curves for n = 4 to 7 can 
be collapsed on to a single curve. The value of x\ obtained by extrapolating 
the curves Zn(x*(e))/n, for n = 7 to 4 is 5.617. We have verified that the 
n = 2 curve and to a lesser extent n = 3 curve deviates from the collapsed 
curve reflecting that nonlinearity corresponding to reinjection is dominant for 
these two cases. This also implies that the orbits do not approach the fixed 
point close enough that the linearized eigenvalues could be useful. Using the 
value of p = 0.0073 and uj = 0.1348 for e = 267.0, we can now estimate the 
value of small period oscillations for n — 11 is 171.0 which is larger than e/2. 
(Here e = 267.0.) This means that a maximum of nine small amplitude oscil- 
lations are allowed beyond the first minimum at this value of e = 267.0. This 
is consistent with what is seen in Fig. 12. We have verified that this method 
of estimating the maximum number of allowed small amplitude oscillations 
for any given value of e works very well as long as n > 3 for m = 1.2 and 
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other values of m as well. For n = 2 and 3, the value of e at which these orbits 
disappear shifts to much lower values that what actually observed. 

A little reflection on the above results reveal the cause of incomplete approach 
to the homoclinic point. Recall that the growth of small amplitude oscillations 
is controlled by 2np/uj. We have seen that this quantity depends inversely on 
e. Thus, the arithmetically increasing number of small amplitude oscillations 
accommodated on Si ( without the trajectory crossing over to S2) is a di- 
rect result of 'softening' of 2irp/uj as a function of e. In other words, for every 
small amplitude oscillation accommodated on Si, e changes by a fixed amount 
commensurate with the softening rate. Thus, the number of small amplitude 
oscillations that can be accommodated in the allowed interval of e in the bifur- 
cation plot will be limited. This also implies that the approach to homoclinic 
point can at best be asymptotic due to finite rate of softening of 2np/u, with 
the asymptotic nature manifesting only in the limit k 2 — > 0. Clearly, these 
results are valid under the assumption that the contribution from nonlinear 
terms to the growth of the small amplitude oscillations is not strong, which is 
substantiated by the numerical evidence that the estimated number of small 
amplitude oscillation allowed for a particular value of e agrees with what is 
numerically observed. 

It may be worth pointing out here that even though the analysis given here 
is for MMOs of the kind I s , it is clear that they can be easily generalized to 
L s kind of MMOs. It must be mentioned here that the 'softening' of 2np/uj is 
a result of the global constraint in our model, namely, the back-to-back Hopf 
bifurcation. Thus, we see that the apparent homoclinic scenario exhibited by 
the model system is completely new. We have shown that this feature coupled 
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with the mechanism of the relaxation oscillations operating in the model gives 
rise to features of MMOs common to both Shilnikov and Gavrilov-Shilnikov 
scenario. 

4-3 Discussion and Conclusions 

Some comments may be in order here about the bent-slow manifold structure 
of the system. The relaxation oscillations seen in this system differs qualita- 
tively from that seen in systems with the S— shaped slow manifold. The major 
difference in the structure is that while in the S— shaped manifold, there are 
two attractive pleats separated by a repulsive part, in our case, both pieces 
Si and 62 of the bent-slow manifold are attractive and are connected continu- 
ously. This aspect coupled with the fact that there is no repulsive part in the 
bent-slow manifold as in the S- shaped manifold suggests that the mechanism 
causing jumps between the Si and S 2 is very different. The only similarity is 
that the number of jumps ( fast transitions) accomplished by the trajectories 
from one part of the slow manifold (£2) to another (Si) and vice versa is two 
as in the case of S— shaped manifold. During the jump from S2 to Si, the 
trajectory tends to move out of the slow manifold into unstable phase space 
by sticking to the direction of of motion on S 2 . Here the motion is acceler- 
ated due to unstable nature of the phase space. This has some similarity to 
canard type of solutions but the comparison is superficial since the unstable 
part of the phase space to which the trajectory moves is not a part of the 
slow manifold. There is another difference namely the operative time scales 
in the dynamics on the slow manifold. The dynamics on S 2 is slow as it is 
controlled by the slow variables y and 0, since x ~ for the entire interval of 
time the trajectory is on S 2 . On the other hand, on Si, the time dependence 
of a trajectory is largely controlled by the fast variable x. 
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We have also analyzed the effect of bent-manifold structure on mixed mode 
oscillations and the incomplete homoclinic scenario. We have shown that in 
the case of bent-manifold structure, the approach of the fast variable towards 
the fixed point is along the slow manifold itself eventhough the eventual ap- 
proach is along the z direction. This is in direct contrast with the S— shaped 
structure where the direction of the jump of the fast variable is transverse 
to the slow manifold containing the fixed point as in the Rossler's S— shaped 
slow manifold. There is another possibility as pointed out by Koper et al\\l 1||. 
In this case, the plane of relaxation oscillations is parallel to the plane of 
nearly harmonic small amplitude oscillations. These authors suggest that this 
type of dynamics may be responsible for the incomplete approach to homo- 
clinicity since the approach to the fixed point is along the attracting pleat 
of the S— shaped manifold. However, in both cases, irrespective of the loca- 
tion of the fixed point, the nature of the relaxation oscillations remain the 
same. Moreover, in both these cases, the small amplitude oscillations(near 
harmonic), as well as the large amplitude (relaxation type) oscillations are 
well characterized. However, the nature of intermediate amplitude oscillations 
are not so well understood. It has been suggested in the literature that these 
are related to canard type of solutions |8|j54|1. The latter type of oscillations 
result from 'sticking' of the trajectory to the repelling part of the S— shaped 
slow manifold before jumping to the attracting pleat of the slow manifold. In 
our case, although the oscillations have a superficial similarity with canard 
type of solutions, it is the 'sticking' of the trajectories in the same direction of 
5*2 well into the unstable phase space, coupled with the fact that there is no 
inherent constraint relating S\ to 5*2 in the manifold structure that appears 
to lead to jumps of all sizes. As an illustration, we have shown a plot of the 
trajectory for m — 1.8 and e = 190.0 in the x — 5 plane (Fig. 14). It is clear 
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that while the small amplitude oscillations are located in the neighborhood of 
the fixed point, the intermediate amplitude oscillations result from the trajec- 
tory 'sticking' to the direction of the S2 plane and moving into the 'unstable' 
part of the phase space by varying amounts each time the trajectory leaves S 2 - 

Although the bent-manifold structure is characteristic of our systems, we be- 
lieve similar structure is likely to be seen in many other models and experi- 
mental systems. Particularly, in chemical kinetics where only binary collisions 
are permitted, models are expected to involve only quadratic or biquadratic 
nonlinearities. Such models are promising candidates to exhibit the bent-slow 
manifold structure. In experimental systems, the dominant signature to look 
for would be the trapping of fast variable at small values over a substantial 
portion of its period followed by sharply peaked pulse like behavior. 

From the above discussion, we see that the bent-slow manifold structure is 
at the root of understanding of the pulsed type of relaxation oscillation and 
the MMOs. We note here that our analysis is completely local since we have 
used the linearized eigenvalues around the fixed point. For the same reason, 
the scaling relation obeyed by the small amplitude oscillations (Eq. 21 and 
Fig. 13) is found to be valid only where the influence of nonlinearity to the 
growth of these small amplitude oscillations (as a function of e) is minimal. 
In spite of this, the scaling relation so obtained forms the basis for estimating 
the maximum number of small period oscillations permitted for a given e thus 
explaining the origin of MMOs in the model. 

Here, we mention that the relaxation oscillations arising out of the atypical 
bent-slow manifold structure is directly related to a dominant characteristic of 
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the PLC effect, namely, the negative strain rate sensitivity of the flow stress. 
The analysis of time scales involved in the relaxation oscillations has been 
useful in understanding the origin of the negative strain rate sensitivity of the 



flow stress which is reported elsewhere [55]. 



In summary, we have analyzed the dynamics of a model for a type of plastic 
instability due to Ananthakrishna and coworkers with particular attention to 
the complex dynamics exhibited by the model. We have shown that the nature 
of relaxation oscillations and the MMO sequences exhibited by the model is 
atypical. We have proposed a new mechanism for the relaxation oscillations 
based on the bent-slow manifold structure of the model. Using this we have 
explained the origin of the MMOs. We have further shown that a crucial 
role in organizing the dynamics is played by the physical constraint, namely, 
the stress oscillations are seen only in a window of strain rates. (Indeed, the 
model has been devised to be consistent with this experimental feature.) This 
constraint translates to back-to-back Hopf bifurcation in the model leading 
to the 'softening' of the eigenvalue of the saddle fixed point that controls the 
small amplitude oscillation. It is this finite rate of softening that is responsible 
for the incomplete approach to the homoclinic bifurcation. 
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Fig. 1. Phase diagram of the model in (m, e) plane. The broken line corresponds to 
the locus of Hopf bifurcations, dotted lines to the PD bifurcations and the continuous 
lines to the locus of SN bifurcations. The thick lines represent the SN bifurcations 
of the PPO culminating in a codimension 2 cusp bifurcation point shown as filled 
diamond. Approximate region of MMOs is shown by the hatched area. 
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Fig. 2. Bifurcation diagram for m = 2.215. 
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Fig. 3. The bifurcation diagram for m = 1.8. Dashed lines indicate unstable periodic 
orbits. Filled circle correspond to a PD bifurcation and filled triangle correspond to 
SN bifurcation. An interior crisis point is marked by an arrow. 
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Fig. 4. Bifurcation diagrams for m = 1.2. The region between the arrows correspond 
to the coexistence region of the MMOs and small amplitude periodic orbits. 
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Fig. 5. Bifurcation diagram for m = 1.0. The dashed line indicates the unstable 
PPO. Inset shows enlarged portion of the bifurcation diagram where the secondary 
Farey sequence manifests immediately after the SN bifurcation of the PPO. 
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Fig. 6. Evolution of a trajectory (thin lines) along with the bent-slow manifold (Si 
and 52 shown by thick lines) structure in the x — 5 plane, for m = 1.2 and e = 200. 
Inset shows the time series of the x variable (continuous line) and z variable (dotted 
line) . 
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Fig. 7. Evolution of the trajectory along with bent-slow manifold (Si and S2) 
structure in (x, y, <p) space indicated by the gray plane, for m = 1.2 and e = 200.0. 
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Fig. 8. Evolution of the trajectory along the bent-slow manifold (Si and S2) struc- 
ture for to = 1.2 and e = 267.0. 
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Fig. 10. Scaling of the period of the superstable orbits with the parameter for 
m = 1.4. The exponential scaling indicates the apparent approach to homoclinicity 
of the The exponential scaling indicates the apparent approach to homoclinicity of 
the saddle focus fixed point. Here e* is taken to be the value of the twelve period 
orbit, e = 247.63. 
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Fig. 11. Plot of a periodic orbit in (x, z, 4>) space for m = 1.2 and e = 267.0 showing 
the eventual direction of approach towards the fixed point is z direction. 
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Fig. 13. Scaling of minimum value of x with e _1 for n = 2, 3, ■ ■ ■ , 7 (m = 1.2). Inset 
shows collapse of all the curves from n = 4 to 7 onto a single curve. Dashed lines 
are shown as guide to the eye. 
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